function vmd_equil
%
%   vmd_potential(fbase, fnum)
%
%   Reads density file. Writes potential and electric field.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Parameters for background potential and system size are edited in the file.


%---------------------------------------------------------

%-- system size

jz=16;
jv=256;

%-- physics paramters

phi00        = 0.005; 
Kz           = 0.34; 
nz           = 1;
Vmax         = 5.0;
vPhase       = 3.5;
ModAmp       = 1.e-9;
Dxxxx        = 1.e-3;
Dvvvv        = 1.e-8;

friction     = 1.e-3;
dt           = 0.1;

%---------------------------------------------------------


%-- set grid vectors --


  dz = 2*pi*nz/(jz*Kz);
  z  = (0:jz-1)*dz;
  dkz = Kz/nz;
  kz = [0:jz/2,-jz/2+1:-1] *dkz;

  dv  = 2*Vmax/jv;
  v   = -vPhase - Vmax + dv*(0:jv-1);
  dkv = 2*pi/(2*Vmax);
  kv  =  [0:jv/2,-jv/2+1:-1] *dkv;

  v0 = -vPhase;

  ExpKz = exp( i * kz * v0 * dt );
 
  ExpKv = zeros(jv,jz);
  for ikz=1:jz
    ExpKv(:,ikz) = exp( -i*v(:)*kz(ikz)*dt - Dxxxx * kz(ikz)**4 *dt);
  end



%-- dealiasing limits

  indz = round(jz/2-jz/6) : round(jz/2+2+jz/6);
  indv = round(jv/2-jv/6) : round(jv/2+2+jv/6);


%-- set background potential and initial conditions

  phi0 = phi00 *exp(i*Kz*z);
  phi0FT = fftn(phi0);


  f0=( 1/sqrt(2*pi) )*exp(-.5*(v + vPhase).**2);
  fvc0T = fft(f0)/jv;


  dfu = zeros(jv,jz);
  for iz=1:jz
    for iv=1:jv
      dfu(iv,iz) = ModAmp * f0(iv) * cos(Kz*z(iz));
    end
  end


 
  for istep=1:1000

  %-- solve poisson equation to find Ez


    rho = sum(dfu,1)*dv;

    phiFT = fft(rho) ./ kz.^2;
    phiFT(1) = 0;

    phiInternal = 2*abs(phiFT(2))/jz;

    phiFT = phiFT + phi0FT;
    EzFT  = -i*kz.*phiFT;  
    Ez  = real(ifftn(EzFT));


    %-- advect in z: df/dt = -v df/dz --


    for iv=1:jv

      fz = dfu(iv,:);

      fz = fft(fz)/jz;
 
      fz(indz)=0;

      fz = fz.*ExpKv(iv,:);

      dfu(iv,:) = ifft(fz);

    end



   %-- acceleration: df/dt = Ez df/dv

   for iz=1:jz

       fvc = dfu(:,iz)';

       fvc = fft(fvc)/jv;

       fvc(indv)=0;

       zz  =  - Dvvvv*kv.**4  - i*kv*Ez(iz) + 1.E-09;

       fvc = fvc.*exp(zz*dt) + (1 - exp(zz*dt)) .* (i*kv*Ez(iz)) .* fvc0T ./ zz;

       dfu(:,iz)= real(ifft(fvc));

   end


   %dfu = dfu * exp( - friction * dt);


   if (rem(istep,100) == 0) 
        printf("%f  %e\n", istep*dt, phiInternal);
   end


end


   fzv =  dfu + meshgrid(f0, 1:jz)';


   %figure(1, "Name", "dfu")
   %imagesc(flipud(dfu));


   %figure(2, "Name", "fzv")
   %imagesc(flipud(fzv))

   dfumax = max(max(abs(dfu)));
   fzvmax = max(max(abs(fzv)));

   one_image(dfu, -dfumax, dfumax, ['dfu.png' ])
   one_image(fzv, -fzvmax, fzvmax, ['fvz.png' ])

end


%---------------------------------------------------------
